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ABSTRACT 

The probability density function (PDF) of the gas density in turbulent supersonic flows is inves- 
tigated with high-resolution numerical simulations. In a systematic study, we compare the density 
statistics of compressible turbulence driven by the usually adopted solenoidal forcing (divergence- free) 
and by compressive forcing (curl- free). Our results are in agreement with studies using solenoidal forc- 
ing. However, compressive forcing yields a significantly broader density distribution with standard 
deviation 3 times larger at the same rms Mach number. The standard deviation-Mach number 
relation used in analytical models of star formation is reviewed and a modification of the existing 
expression is proposed, which takes into account the ratio of solenoidal and compressive modes of the 
turbulence forcing. 



Subject headings: hydrodynamics — ISM: clouds — ISM: kinematics and dynamics ■ 
— methods: numerical — turbulence 
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1. INTRODUCTION 



The pioneering works by [P adoan et al.' (1997", PNJ97 
below) and iPassot fc Vazqu ez-Seniadeni (1998, PV98 
below) have shown that the standard deviation (stddev) 
(Tp, i.e., the width or dispersion of the linear density 
PDF pp grows proportional to the rms Mach number 
M of the turbulent flow, however, with proportionality 
constants in disagreement by a factor of 2. The exact 
dependence of the PDF's stddev on the rms Mach num- 
ber is a key ingredient for analytical models of star for- 
mation. For instance, [Padoan fc Nordlundl (2002) and 
iHennebelle fc Chabrieil ()2008D relate the densitv PDF to 
the core mass f unction (CMF) and stellar initial mass 
function (IMF). |Tassii (j2007t ) uses the density PDF on 
galac tic scales to repro duce the Kennicutt-Schmidt rela- 
tion. lElmegreenl (|2008D suggests that the star formation 
efficiency is a function of the density PDF. 

In the present work, we solve the discrepancy between 
PNJ97 and PV98 by showing that the stddev of the PDF 
is not only a function of the rms Mach number, but is 
also very sensitive to the relative importance of solenoidal 
(divergence-free) and compressive (curl-free) modes of 
the turbulence forcing, leading to variations in the std- 
dev up to factors of ~3 for the same rms Mach number. 
The main result of the present work is that the conclu- 
sions of PNJ97 and PV98 can be harmonized, if we take 
into account that PNJ97 have analyzed purely solenoidal 
forcing, whereas PV98 have studied purely compressive 
forcing. This apparent difference has not been considered 
analytically or numerically before. 

We begin by explaining our numerical method in ij2l ij3] 
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shows that our results are consistent with previous stud- 
ies using solenoidal forcing, whereas compressive forc- 
ing yields a PDF with stddev 3 times larger in 3- 
dimensional turbulent flows. We present a heuristic 
model explaining this difference, which is based on the 
ratio of solenoidal to compressive modes of the forcing 
to estimate the proportionality constant in the stddev- 
Mach relation. 21 summarizes our conclusions. 

2. NUMERICAL METHODS 



Th e piecewise parabolic method (jColella fc WoodwardI 
1984") implementation of FLASH3 (iFryxell et all l2000l : 
Dubey et al. 2008) was used to solve the hydrodynamic 
equations on periodic uniform grids with 16384^, 4096^, 
and 1024^ grid points in 1, 2 and 3 dimension(s) (ID, 2D, 
3D). Since we model isothermal gas, it is convenient to 
define s = ln(p/po) as the natural logarithm of the den- 
sity divided by the mean density po- Isothermal gas has 
been modelled several times in the context of molecular 
cloud dynamics taking periodic boundaries and studying 
compressible turbulen ce with solenoidal or weakly com- 
pressivc forcing (e.g.. Stone et al. 1998; Mac Low et al 
1998: Mac Low 1999- Klcsscn ct al. 2000: Hcitsch et al 
2001[ iBoldvrev etHl 120021: ILi et all 12 003: Padoan et al 
2QO4I; 'Ja pDsen et all 120051: iBallest eros-Parcdes et al 
2006; Krit suk et aLn2007l : iDib et all 2008; Offner et al 
2008). Most of these studies have purely solenoidal 
forcing motivated by incompressible turbulence studies, 
some have weakly compressive forcing. The latter is a 
result of the natural ratio of solenoidal to compressible 
modes of a Gaussian field prepared in Fourier space with- 
out subsequent projection. This natural ratio is 2:1 in 3D 
and 1:1 in 2D. Accordingly, the ratio of compressive to 
the sum of compressive plus solenoidal modes is 1:3 in 
3D and 1:2 in 2D. Only PV98 have purely compressive 
forcing, because they analyze ID simulations where no 
solenoidal component exists and the ratio of compressive 
to total modes is 1:1. 

The forcing is typically eithe r modelled as a static 
patt ern following the r ecipes by iMac Low et al.l (|1998l ) 
and IStone et al.l (|1998l ). or by an Ornstein-Uhlenbeck 
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Fig. 1. — Column density field in units of the mean column density for solenoidal forcing (left), and compressive forcing (right) at a 
randomly picked time in the regime of statistically stationary turbulence. Both maps show 4 orders of magnitude in column density with 
the same scaling for direct comparison of solenoidal and compressive forcing at rms Mach number ~5.5. 



(OU) process with fin ite autocorrelation timescale T 
(e.g., EswaranXPoEd[l988; Schmidt ct al. 2006). Note 
that driving with a static pattern is the hmiting case of 
an OU process with infinite autocorrelation timescale. 
We use the OU process and follow the usual approach, 
setting T equal to the dynamical timescale T = L/{2V) 
{L is the size of the computational domain, V = CsM. 
and is the rms Mach number) , which is equal to 
the decay time const ant of the turbulence (jStone et al.1 
Il998t iMac Lowl Il999l ) at the scales of energy injection 
1 < A; < 3. This guarantees a well-defined stochastic 
driving field, that varies smoothly in space and time. 
We checked that modelling of the forcing as an almost 
static pattern (increasing T by one order of magnitude) 
did not significantly affect our results. 

In order to obtain a purely solenoidal, or a purely 
compressive forcing (or any combination), a Helmholtz 
decomposition can be made by applying the projection 
operator in Fourier space (fc-space) 

n^- = cn^ + (1 - c)4 = c% + (1 - 2C)|^ , (1) 

to the random vector field returned by the OU process 
or generated by the usual recipes. The parameter Q G 
[0, 1] controUs the relative importance of solenoidal and 
compressive modes. If we set ^ = 1, V'lj projects only 
the solenoidal component, whereas only the compressive 
component is obtained by setting C = 0. 

3. RESULTS AND DISCUSSION OF THE STANDARD 
DEVIATION - MACH NUMBER RELATION 

Figure [T] compares column density maps for solenoidal 
vs. compressive forcing of our 3D models (1024'^) from 
a randomly picked snapshot in the regime of statisti- 
cally stationary turbulence. Obviously, compressive forc- 
ing yields much larger density contrasts, despite the fact 
that the rms Mach number ^5.5 is about the same in 
both cases. This is quantitatively shown in Fig. [21 which 
presents the comparison of the time averaged volume- 
weighted density PDFs in the solenoidally and compres- 




FlG. 2. — Volume-weighted density PDFs p(s) in linear scaling 
where s = ln(p/po). The PDF obtained by compressive forcing 
(comp, f = 0.0) is much broader compared to the solenoidal one 
(sol, C, = 1.0) at the same rms Ma ch numb er. The p eak is shifted 
due to mass conservation l lVazau ez-Scmad enil 11994 ). Gray error 
bars indicate 1-sigma temporal fluctuations of the PDF. A sample 
of ~10^^ datapoints contribute to each PDF. 

sively driven cases. Most importantly, although the rms 
Mach numbers are almost the same, the stddev Cp of 
the PDF obtained in compressive forcing is ~ 3 times 
larger than for solenoidal forcing (Tab. [T|). This is the 
main result of the present study and has important con- 
sequences: The stddev of the density PDF is not only 
a function of the rms Mach number as found by PNJ97 
and PV98, but also a function of the relative strength of 
solenoidal to compressive modes of the turbulence forc- 
ing, i.e., a function of the projection parameter C, in 
eq. ([T]). In the following, we review the stddev-Mach 
number relation and present a heuristic model for the 
proportionality constant in this relation. 

It has been pointed out by PNJ97 and measured by 
PV98 that the stddev ap of the PDF of the linear density 
Pp grows linear with the rms Mach number M. like 

(Jp/pa = bM (2) 
with the proportionality constant b. PV98 find 6 ^ 1 in 
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TABLE 1 

Parameters in solenoidal and compressive forcing 



D 


resolution 


c 


M 


O"p/P0 


6 


3 


1024^ 


1.0 (sol) 


5.3±0-2 


1.9±0-i 


0.36±0-03 


3 


1024^ 


0.0 (comp) 




5.9±io 


1.05±0-i9 


2 


4096^ 


1.0 (sol) 


5g±0.5 


2,4±0.5 


0.43±o.io 


2 


4096^ 


0.0 (comp) 


5 7±0-6 


5.0±i-5 


0.88±0-28 


1 


163841 


0.0 (comp) 


5.0±io 


4.5±2.o 


0.90±0-44 



ID simulations with M ranging from subsonic to super- 
sonic flows. PNJ97 motivate and explain eq. ^ with the 
isothermal hydrodynamic shock jump conditions. They 
assume the PDF follows a lognormal analytical form 



Ps{s) 



exp 



2(ji 



(3) 



to get an expression for the stddev in the logarithmic 
density 

[ln(l + 6W)]'/' . (4) 

Using Psds ~ Ppdp, it is easy to show that for 
any density PDF, whether lognormal or not, cr^ = 

Pq J^^[ex.p{s) — 1]^ Psds. From the assumption ^ it 
follows that eq. ([4]) and ^ are equivalent and that both 
expressions have the same b. This means that eq. ([2]) 
is the basic stddev-Mach number relation from which 
eq. (|4]) follows, if a lognormal PDF is assumed. 

PNJ97 applied eq. ^ to magnetohydrodynamical 
(MHD) simulations and obtained &~0.5. Consequently, 
their stddev i s only h alf as large as that found by PV98. 
iPadoan etU (2007*) argue that eq. ^ may hold for 
MHD turbulence as well, if A4 is replaced by the Alfvenic 
Mach number, although the shock jump conditions mo- 
tivating eq. ^ are different for shocks perpe ndicular 

to the mag netic fiel d. In MHD simulations, iLi et al.l 

(|200l and iLi et all ([2008^ find b - 0.38 and b ~ 0.41 
for ideal MHD and &~0.58 by including ambipolar dif- 
fusion (AD). The latter means that including AD broad- 
ens the PDF over ideal MHD, because of the reduced 
magnetic pressure compared to ideal MHD. Therefore, 
the parameter b in eq. ([2]) and ^ is not expected to 
be universal, b ut dep ends on the magnetic field (see 
also.lVazauez-S emadcn r"et al.ll2005l) . On the other hand, 
iLemaster fc~Stona (.2008^ find that even strong magnetic 
fields do not alter the stddev-Mach number relation sig- 
nificantly. In the following discussion, we will concen- 
trate on purely hydrodynamical estimates showing that 
b is much more sensitive to the way of forcing compared 
to the reported MHD effects. In order to show this, we 
use the more fundamental stddev-Mach number relation 
given by eq. ([2]), which does not rest on the additional 
assumption of a lognormal density PDF. 

From purely hydrodynamic studies utilizing different 
numerical methods and resolutions, typically smaller but 
roughly consistent values of b are found compared to 
6 ~ 1 by PV98, e.g^^, b ~ 0.35 using 200,000 SPH par- 
ticles (Li et al. 200lj_6~ 0.26 using ENZO with 1024^ 
cells (|Kritsuk et alTf2007D '^. b 0.37 using anothe r grid- 
based method with 512^ cells ()Beetz et al.ll2008| ). It is 

* In their recent 2048'^ study with solenoidal forcing, they find 
a slightly larger 6~0.32 (Kritsuk 2008, private communication). 
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Fig. 3. — Density PDFs p(s) obtained from 3D, 2D and ID sim- 
ulations with compressive forcing and from 3D and 2D simulations 
using solenoidal forcing in logarithmic scaling. Note that in ID, 
only compressive forcing is possible as in the study by PV98. As 
suggested by eq. l(5)l, compressive forcing yields almost identical 
density PDFs in ID, 2D and 3D with 6 ~ 1 because compressive 
forcing induces gas compression along all available spatial direc- 
tions. Solenoidal forcing on the other hand leads to a density PDF 
with fe~l/2 in 2D and with b'^1/3 in 3D as a consequence of the 
natural ratio of compressible to compressible plus solenoidal modes 
building up in the velocity field, which is 1:2 in 2D and 1:3 in 3D 
for solenoidally driven supersonic turbulent flows. 

important to note that all the aforementioned studies (in- 
cluding the MHD studies) use solenoidal or weakly com- 
pressive forcing, except for PV98, who naturally have 
purely compressive forcing because they analyze ID sim- 
ulations. Similar to PV98, we find b ~ 1.05 for purely 
compressive forcing (Tab.[T]). For purely solenoidal forc- 
in g, we get & ^ 0. 36 in agreement wi th t he estimates 
by Li et all (|2003D . lKritsuk et al.l (|2007D and lBeetz et~aLl 
(200a ). This shows that the width of the PDF does not 
only depend on the rms Mach number given by eq. ([2]), 
but also on the mixture of solenoidal and compressive 
modes ( of the forcing. This dependence was q uali- 
tatively mentioned by Nordhmd & Padoan (1999), al- 
though up to now, no quantitative estimate existed. 

Based on the diversity of b in eq. dH) and dH) obtained 
in the studies mentioned above and based on our di- 
rect comparison of solenoidal and compressive forcing, 
we suggest that a plausible estimate for the proportion- 
ality constant b is obtained by taking into account the 
compressibility induced by the forcing. 



b = 1 



1-1 
D 



c = 



1 - (2/3)C, 
1-(1/2)C, 
1, 



for 
for 
for 



D = 3 
D = 2 



(5) 



The relative strength of solenoidal to compressive modes 
C is defined in the projection operator, eq. ([1]). 

For solenoidal forcing (C = 1), by construction no com- 
pression is directly induced by the forcing. Rather, den- 
sity fluctuations are solely a result of compressions nat- 
urally building up in a supersonic turbulent flow. This 
is a consequence of the ratio of compressible modes to 
the sum of compressible plus solenoidal modes in the ve- 
locit y field, that the turbulent gas adjusts itself to (see, 
e.g.. lElmegreen &: Scald [20041 ) . This ratio is related to 
the number of degrees of freedom D (spatial directions) 
available for compressible modes and is simply given by 
1/D showing up in eq. Consequently, eq. (O leads to 
6^^1/3 for the 3D case with solenoidal forcing. This is in 
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reasonable agreemen t with the res u lts bylLi et aTl (|2003D , 
iKritsuk et all ()2007f ). iBeetz et alJ (|2008h and ours. 

The ID case as analyzed by PV98 is a special case of 
compressive forcing in arbitrary dimensions. Note that 
for compressive forcing, eq. ^ always suggests 6 ~ 1 
independent of D. This can be understood as a conse- 
quence of the direct compression induced by the compres- 
sive force field, contrary to solenoidal forcing for which 
the natural ratio of modes in the velocity field determines 
b. Compressive forcing by construction immediately in- 
duces compressions along all available spatial directions 
D directly, and b ^ D/D = 1. The ID calculations 
(C = 0) by PV98 therefore exhibit ~ 1. Similar, in 
our high-resolution 3D case of purely compressive forcing 
(C = 0), compression is induced along all the 3 available 
spatial directions resulting in 6~3/3 = 1. 

Besides testing eq. ([SJ for the extreme cases of 
solenoidal and compressive forcing at different dimen- 
sionality, we also aim at testing it for intermediate val- 
ues of C (mixtures of solenoidal and compressive forcing) 
using reported results in the literature and additional 
n umerical simulation s. 

I Vaz quez- Semade"nil (|1994l ) have analyzed 2D simula- 
tions with a 1:1 mixture of solenoidal and compres- 
sive forcing — 0.5). Thus, the expected value is 
6 ~ 0.75, and they fin d b^O.7 close to the expectation. 
I Schmidt et al.l (|2008l ) use a mainly compressive forcing 
(C = 0.1) obtaining 6~0.94 from their data, while eq. ([5]) 
suggests 6~0.93. 

We ran additional 2D and ID simulations with 
solenoidal and compressive forcing to support our heuris- 
tic model, eq. ([5]). In order to get a statistically signif- 
icant sample, the numerical resolution along the spatial 
directions was increased from 1024 in 3D to 4096 in 2D 
and 16384 in ID, as well as the integration times were 
increased yielding at least twice as many sampling points 
as PV98 have. Table [1] summarizes the parameters of all 
simulations and numerical estimates for b. As expected 
from eq. (O, the compressive cases in ID, 2D and 3D, 
all exhibit nearly the same & ~ 1. Moreover, the PDFs 
are very similar in all of the compressively driven cases, 
which is shown in Fig. [3l The 3D case with solenoidal 
forcing is also included for comparison (6~l/3), as well 
as the 2D case with solenoidal forcing, for which we es- 
timate 6^^0.43 close to the prediction 6^1/2. 

4. SUMMARY AND CONCLUSIONS 

Performing high-resolution hydrodynamical simula- 
tions of driven isothermal compressible turbulence, we 
have found that the way of forcing the turbulence has 
a strong effect on density statistics. We have compared 
the usually adopted solenoidal forcing (divergence-free) 
and compressive forcing (curl- free). The most important 
result is that in 3D, compressive forcing yields a density 
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PDF with stddev ^ 3 times larger compared to solenoidal 
forcing for the same rms Mach number. 

As found by PNJ97 and PV98, the stddev cTp is increas- 
ing directly proportional to the rms Mach number, how- 
ever, they find different proportionality constants. Tak- 
ing into account the ratio of solenoidal to compressive 
modes C of the forcing resolves the disagreement between 
PNJ97 and PV98. We suggest that the proportional- 
ity constant b in the stddev-Mach number relations ([2]) 
and ^ can be determined by a heuristic model sum- 
marized in eq. ([5]). In the case of compressive forcing 
(C = 0), the proportionality constant &~1 irrespective of 
the dimensionality of the simulation. Solenoidal forcing 
on the other hand yields 6'--^ 1/3 in 3D and 6^1/2 in 2D. 

We mention the impact of our results on analyti- 
cal models linking the stat istics of supersonic tu rbu- 
lence to the CMF/I MF. Pad oan fc Nordlundl (|2002D and 
iHennebelle fc Chabrier. (,2008f) rely on integrals over the 
density PDF to get a handle on the mass of objects 
above a certain density threshold. Since the width of 
the PDF is so sensitive to the mixture of solenoidal 
and compressive modes Ci we also expect a strong in- 
fluence on the derived CMF/IMF. Indeed, the larger 
dispersion obtained from compressive forcing leads to 
better agreement (Hennebelle 2008, private communica- 
ti on) of the analytic expr ession with the observed IMF 
in IHennebelle fc Chabrie^ (j2008). 

Given that many pro posed sources of inter- 
stellar turbulence (e.g., lElmegreen fc Scalol 120041 : 
iMac Low fc KlessenI l2004f ) are likely to directly ex- 
cite compressive modes (e.g., galactic spiral density 
shocks, large scale gravitational contraction, supernova 
explosions, protostellar jets, winds and outflows), it 
is reasonable to expect that turbulence in the ISM 
is driven by a mixture of solenoidal and compressive 
modes, possibly with compressive modes being more 
important than solenoidal modes. 
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